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ABSTRACT 

We present a new analysis of the ionizing emissivity (N- Ion , s -1 Mpc -3 ) for galaxies dur¬ 
ing the epoch of reionization and their potential for completing and maintaining reion¬ 
ization. We use extensive SED modelling incorporating two plausible mechanisms 
for the escape of Lyman continuum photon - to explore the range and evolution of ion¬ 
izing efficiencies consistent with new results on galaxy colours (/3) during this epoch. 
We estimate N[ on for the latest observations of the luminosity and star-formation rate 
density at 2 < 10, outlining the range of emissivity histories consistent with our new 
model. Given the growing observational evidence for a UV colour-magnitude relation 
in high-redshift galaxies, we find that for any plausible evolution in galaxy properties, 
red (brighter) galaxies are less efficient at producing ionizing photons than their blue 
(fainter) counterparts. The assumption of a redshift and luminosity evolution in (3 
leads to two important conclusions. Firstly, the ionizing efficiency of galaxies natu¬ 
rally increases with redshift. Secondly, for a luminosity dependent ionizing efficiency, 
we find that galaxies down to a rest-frame magnitude of Muv ~ —15 alone can poten¬ 
tially produce sufficient numbers of ionizing photons to maintain reionization as early 
as 2 ~ 8 for a clumping factor of Chu < 3. 

Key words: dark ages, reionization, first stars - galaxies: high-redshift - galaxies: 
luminosity function, mass function - galaxies: evolution 


1 INTRODUCTION 

At the present day, the intergalactic and interstellar medium 
(IGM,ISM) are known to be predominantly ionized. How¬ 
ever, following recombination at z fa 1100, the baryon con¬ 
tent of the Universe was mostly neutral. At some point in 
the history of the Universe, the IGM underwent a transition 
from this neutral phase to the ionized medium we see today, 
a period known as the epoch of reionization (EoR hereafter). 
The strongest constraints on when reionization occurred are 
set by observations of the Gunn-Peterson trough of distant 
quasars (Fan et al. 2006), which indicate that by 2 ~ 5.5, the 
Universe was mostly ionized (with neutral fractions ~ 10 -4 ). 

Additionally, measurements of the total optical depth 
of electrons to the surface of last scattering implies that 
reionization should be occurring at higher redshift, towards 
z « 10, for models of instantaneous reionization (Hinshaw 
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et al. 2013; Bennett et al. 2013). However, critical outstand¬ 
ing questions still remain. Firstly, when did the intergalactic 
hydrogen and helium complete reionization? And secondly, 
what were the sources of ionizing photons which powered 
the reionization process? Was it predominantly powered by 
star-forming galaxies or by active galactic nuclei/quasars? 

For hydrogen reionization with its earlier completion, 
the rapid decline in the quasar luminosity function at high 
redshift (Willott et al. 2010; Fontanot et al. 2012, 2014) 
does suggest that star-forming galaxies are the most likely 
candidates for completing the bulk of reionization by 2 ~ 
6. The contribution from faint AGN could still however 
make a significant contribution to the ionizing emissivity 
at 2 > 4 (Giallongo et al. 2015). Based on the optical depth 
constraints set by WMAP (Hinshaw et al. 2013) and ei¬ 
ther observed IGM emissivities at lower redshift (Kuhlen 
& Faucher-Giguere 2012; Robertson et al. 2013; Becker & 
Bolton 2013), or emissivities predicted by simulations (Cia- 
rdi et al. 2012), several studies have drawn the same con- 
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elusion that faint galaxies from below the current detection 
limits and/or an increasing ionizing ‘efficiency’ at higher red- 
shift is required. Even with the lower optical depth measure¬ 
ment now favoured by the recent Planck analysis (Planck 
Collaboration et al. 2015), such assumptions are essentially 
still required to satisfy these criteria (Robertson et al. 2015). 

One of the possible mechanisms for this increasing ion¬ 
izing efficiency is an evolution in the fraction of the ionizing 
photons able to escape their host galaxy and ionize the sur¬ 
rounding IGM, known as the escape fraction (/ esc ). There 
have been several studies designed to understand this issue, 
but there are still large uncertainties in what the escape 
fraction for galaxies is and how it evolves with redshift and 
other galaxy properties. In a study of z ~ 1.3 galaxies, Siana 
et al. (2010) searched for Lyman-continuum photons from 
star forming galaxies, although no systems were detected. 
After correcting for the Lyman-break and IGM attenuation 
the limit placed on the escape fraction is f esc < 0.02 after 
stacking all sources. Bridge et al. (2010) find an even lower 
limit of /esc < 0.01 using slitless spectroscopy at z ~ 0.7, 
although one AGN in their sample is detected. However, 
higher escape fractions of ~ 5% to ~ 20 — 30% have been 
measured for galaxies at 2 ~ 3 (Shapley et al. 2006; Iwata 
et al. 2009; Vanzella et al. 2010; Nestor et al. 2013), consis¬ 
tent with the relatively high / eac ~ 0.2 expected from IGM 
recombination rates determined from Lya forests (Bolton 
& Haehnelt 2007). Furthermore, the average / eB c for galax¬ 
ies at 2 ~ 3 may be significantly higher than the existing 
measurements due to the selection biases introduced by the 
Lyman-break technique (Cooke et al. 2014). 

It is important to bear in mind that the property which 
is fundamental to studies of reionization is the total number 
of ionizing photons which are available to ionize the inter- 
galactic medium surrounding galaxies. Hence, while the es¬ 
cape fraction is a critical parameter for reionization, it must 
be measured or constrained in conjunction with the under¬ 
lying continuum emission to which it applies. For example, 
an increase in / esc may not have an effect on reionization if 
it is accompanied by a reduction in the intrinsic number of 
ionizing photons being produced. 

As shown in Robertson et al. (2013) (see also Leitherer 
et al. 1999), the number of ionizing photons produced per 
unit UV luminosity emitted (e.g. L 1500j O can vary signif¬ 
icantly as a function of the stellar population parameters 
such as age, metallicity and dust extinction. Thankfully, evo¬ 
lution or variation among the galaxy population in these pa¬ 
rameters will not only influence the production of ionizing 
photons but will have an effect on other observable proper¬ 
ties such as the UV continuum slope (/3, Calzetti et al. 1994). 
With the advent of ultra-deep near-infrared imaging surveys 
such as the UDF12 (Koekemoer et al. 2013) and CANDELS 
(Grogin et al. 2011; Koekemoer et al. 2011) surveys, obser¬ 
vations of the UV continuum slope extending deep into the 
epoch of reionization are now available. Furthermore, there 
is now strong evidence for an evolution in /3 as a function of 
both galaxy luminosity and redshift out to 2 ~ 8 (Bouwens 
et al. 2014b; Rogers et al. 2014). 


In this paper we use the latest observations of /3 span¬ 
ning the EoR combined with SED modelling, incorporat¬ 
ing the physically motivated escape mechanisms to explore 
what constraints on the key emissivity coefficients are cur¬ 
rently available. We also explore the consequences these con¬ 
straints may have on the EoR for current observations of the 
star-formation rates in this epoch. In addition to the obser¬ 
vations of in-situ star-formation through the UV luminosity 
functions, we also investigate whether recent measurements 
of the galaxy stellar mass function and stellar mass density 
at high redshift (Duncan et al. 2014; Grazian et al. 2014) can 
provide additional useful constraints on the star-formation 
rates during EoR (Stark et al. 2007; Gonzalez et al. 2010). 

In Section 2, we outline the physics and critical param¬ 
eters required to link the evolution of the neutral hydrogen 
fraction to the production of ionizing photons. We also ex¬ 
plore plausible physical mechanisms for the escape of Lyman 
continuum photons from galaxies, outlining the models ex¬ 
plored throughout the paper. We then review the current 
literature constraints on the UV continuum slope, /3, both 
as a function of redshift and galaxy luminosity. Next, in Sec¬ 
tion 3, we explore how the escape fraction, dust extinction 
and other stellar population parameters affect the observed 
f3 and the coefficients relating ionizing photon production 
rates to observed star-formation rates and UV luminosities. 
In Section 4, we apply these coefficients to a range of exist¬ 
ing observations, exploring the predicted ionizing emissivity 
throughout the epoch of reionization for both constant and 
redshift dependent conversions. We then discuss how the 
varying assumptions and proposed relations would impact 
the ionizing photon budget consistent with current obser¬ 
vations before outlining the future prospects for improving 
these constraints in Section 5. Finally, we summarise our 
findings and conclusions in Section 6. 

Throughout this paper, all magnitudes are quoted in 
the AB system (Oke & Gunn 1983). We also assume a A- 
CDM cosmology with Ho = 70 kms _1 Mpc _1 , = 0.3 
and fi a = 0.7. Quoted observables (e.g. luminosity density) 
are expressed as actual values assuming this cosmology. We 
note that luminosity and luminosity based properties such as 
stellar masses and star-formation rates scales as h~ 2 , whilst 
densities scale as h 3 . 


2 LINKING REIONIZATION WITH 
OBSERVATIONS 

2.1 The ionizing emissivity 


The currently accepted theoretical picture of the epoch of 
reionization, as initially described by Madau et al. (1999), 
outlines the competing physical processes of ionization of 
neutral hydrogen by Lyman continuum photons and recom¬ 
bination of free electrons and protons. The transition from 
a neutral Universe to a fully ionized one can be described 
by the differential equation: 


Qhii 


Vic 


(n H > 


Qhii 

(^rec) 


(1) 
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where Qhh is the dimensionless filling factor of ionized hy¬ 
drogen (such that Qhh = 1 for a completely ionized Uni¬ 
verse) and N \ on is the comoving ionizing photon production 
rate (s -1 Mpc -3 ) or ionizing emissivity. The comoving den¬ 
sity of hydrogen atoms, (tih), and average recombination 
time (tree) are redshift dependent and is dependent on the 
primordial Helium abundance, IGM temperature and cru¬ 
cially the inhomogeneity of the IGM, parametrised as the 
so-called clumping factor Chii = (nl) / (ne) 2 (Pawlik et al. 
2009). We refer the reader to Madau et al. (1999), Kuhlen 
& Faucher-Giguere (2012) and Robertson et al. (2010, 2013) 
for full details on these parameters and the assumptions as¬ 
sociated. 

In this paper we will concentrate on Mon and the pro¬ 
duction of Lyman continuum photons by star-forming galax¬ 
ies. The link between the observable properties of galaxies 
and the ionizing photon rate, Mi on can be parametrised as 

Mion — /esc^ion PU V (2) 

following the notation of Robertson et al. (2013) (R13 here¬ 
after), where puv is the observed UV (1500A) luminosity 
density (in erg s -1 Hz -1 Mpc -3 ), £i on the number of ion¬ 
izing photons produced per unit UV luminosity (erg -1 Hz) 
and / esc the fraction of those photons which escape a host 
galaxy into the surrounding IGM. Alternatively, -M on can be 
considered in terms of the star-formation rate density, psfr 
(M 0 yr -1 Mpc -3 ), 

Vjori • /escKionPSFR (3) 

where Ki on (s -1 Mg 1 yr) is the ionizing photon production- 
rate per unit star-formation (£q in the notation used in 
Robertson et al. (2010)). In addition to observing puv or 
Psfr during the EoR, accurately knowing f esc , £ion and Ki on 
is therefore critical to estimating the total ionizing emissiv¬ 
ity independent of how well we are able to measure the UV 
luminosity or star-format ion rate densities. 

As the production of LyC photons is dominated by 
young UV-bright stars, the rate of ionizing photons is there¬ 
fore highly dependent on the age of the underlying stellar 
population and the recent star-formation within the galaxy 
population. Physically motivated values of £i on (or its equiv¬ 
alent coefficient in other notation) can be estimated from 
stellar population models based on plausible assumptions on 
the properties of high-redshift galaxies (Bolton & Haehnelt 
2007; Ouchi et al. 2009; Kuhlen & Faucher-Giguere 2012). 

However, with observations of the UV continuum slope, 
f3 (Calzetti et al. 1994), now extending deep into the epoch 
of reionization, limited spectral information is now available 
for a large sample of galaxies. Despite the many degeneracies 
in /3 (see later discussion in Section 3), it is now possible to 
place some constraints on whether the assumptions made are 
plausible. In R13, values of £i on are explored for a range of 
stellar population parameters relative to the /3 observations 
of Dunlop et al. (2013). Based on the range of values con¬ 
sistent with the observed values of /3 ~ —2, Robertson et al. 
choose a physically motivated value of log 10 £i on = 25.2. 

Typically, a constant / esc is applied to all galaxies in 


addition to the estimated or assumed values of £ion/Ki on 
(Ouchi et al. 2009; Finkelstein et al. 2012b; Robertson et al. 
2013), motivated in part by our lack of understanding of 
the redshift or halo mass dependence of / eac ■ However, ap¬ 
plying a constant f esc does not take into account exactly 
how the Lyman continuum photons escape the galaxy, what 
effect the different escape mechanism might have on the ob¬ 
served galaxy colours, and how that might alter the assumed 
£ion/ftion based on /3. 

2.2 Mechanisms for Lyman continuum escape 

In Zackrisson et al. (2013), detailed SSP and photo¬ 
dissociation models were used to explore how future ob¬ 
servations of /3 and Ha equivalent-width with the James 
Webb Space Telescope (JWST) could be used to constrain 
the escape fraction for two different Lyman continuum es¬ 
cape mechanisms (see Fig. 8 of Zackrisson et al. (2013)). 
However, it may already be possible to rule out significant 
parts of the / esc and dust extinction parameter space using 
current constraints on /3 and other galaxy properties. To es¬ 
timate the existing constraints on / esc , Cion, Kion and their 
respective products, we combine the approaches of Zackris¬ 
son et al. (2013) and Robertson et al. (2010, 2013). To do 
this, we model /3, Cion and Ki on as a function of f esc for the 
two models of Zackrisson et al. (2013). The components and 
geometry of these two models are illustrated in Fig. 1. 

In the first model, model A hereafter (Fig. 1 left) and 
dubbed ‘ionization bounded nebula with holes’ by Zackris¬ 
son et al., Lyman continuum photons along with unattenu¬ 
ated starlight are able to escape through low-density holes 
in the neutral ISM. In this model, the escape fraction is de¬ 
termined by the total covering fraction of the neutral ISM. 
Under the assumption that these holes are small and evenly 
distributed, the observed galaxy SED (averaged across the 
galaxy as in the case of photometry) would then be a com¬ 
bination of the unattenuated starlight from holes and the 
dust reddened starlight and nebular emission from the Hi 
enshrouded regions. 

A second model, model B hereafter (Fig. 1 right) corre¬ 
sponds to the ‘density bounded nebula’ of Zackrisson et al. 
(2013). This model could occur when the local supply of Hi 
is exhausted before a complete Stromgren sphere can form, 
allowing Lyman continuum photons to escape into the sur¬ 
rounding ISM. The fraction of LyC photons which can es¬ 
cape the nebular region is determined by the fraction of the 
full Stromgren radius at which the nebular region is trun¬ 
cated. The total escape fraction is then also dependent on 
the optical depth of the surrounding dust screen. Of these 
two mechanisms, the former (Model A: ionization bounded 
nebula with holes) is the model which most closely repre¬ 
sents the physics predicted by full radiation hydrodynamical 
models of dwarf galaxies. 

In Wise & Cen (2009), it was found that Lyman contin¬ 
uum radiation preferentially escaped through channels with 
low column densities, produced by radiative feedback from 
massive stars. The resulting distribution of LyC escape frac- 
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Figure 1. Schematic cartoon illustrations of the Lyman continuum escape mechanisms outlined in Section 2. For both models, the stars 
represent a central galaxy surrounded by a Hu ionization region, dust is distributed in an outer dust-screen. Left: An ionization bounded 
nebula with holes (sometimes referred to as the ’picket fence model’) in which LyC escapes through holes in the ISM. Right: a density 
bounded nebula where LyC is able to escape due to the incomplete Stromgren sphere formed when the galaxy depletes its supply of 
neutral hydrogen. 


tion within a galaxy is highly anisotropic and varies signif¬ 
icantly between different orientations. Evidence for such an 
anisotropic escape mechanism has also been found recently 
by Zastrow et al. (2013), who find optically thin ionization 
cones through which LyC can escape in nearby dwarf star- 
bursts. Similarly, Borthakur et al. (2014) find a potential 
high-redshift galaxy analog at 2 ~ 0.2 with evidence for 
LyC leakage through holes in the surrounding neutral gas 
with an escape fraction as high as f esc ~ 0.2 (21%). 

However, this value represents the optimum case in 
which there is no dust in or around the low-density chan¬ 
nels (corresponding to Model A). For the same system, when 
Borthakur et al. (2014) include dust in the low-density chan¬ 
nels, the corresponding total LyC escape fraction is reduced 
to « 1%. The two models explored in this work represent 
the two extremes of how dust extinction will effect the es¬ 
caping Lyman continuum for toy models such as these, the 
dust-included estimates of Borthakur et al. (2014) therefore 
represent a system which lies somewhere between Models A 
and B. 

A potential third mechanism for Lyman continuum es¬ 
cape was posited by Conroy & Kratter (2012), whereby ‘run¬ 
away’ OB stars which have traveled outside the galaxy cen¬ 
tre can contribute a significant amount to the LyC emit¬ 
ted into the surrounding IGM. For high-redshift galaxies 
with significantly smaller radii than local galaxies, massive 
stars with large velocities could venture up to 1 kpc away 
from their initial origin into regions with low column den¬ 
sity. Conroy & Kratter (2012) estimate that these stars could 
in fact contribute 50 — 90% of the escaping ionizing radia¬ 
tion. In contrast, recent work by Kimm & Cen (2014) finds 
that when runaway stars are included into their models of 
Lyman continuum escape, the time average escape fraction 
only increases by ~ 20%. Given the additional complica¬ 
tions in modelling the relevant observational properties and 
their relatively small effect, we neglect the contribution of 
runaway stars in the subsequent analysis. 


In Section 3, we describe how we model both the ob¬ 
servable (/?) and unobservable (£ion,Kaon) properties for both 
model A and model B. But first, we examine the existing 
observations on the evolution of /3 into the epoch of reion¬ 
ization. 

2.3 Observed UV Continuum Slopes 

In Fig. 2, we show a compilation of recent results in the 
literature on the observed UV slope, /?, as a function of 
both redshift and rest-frame UV magnitude, Muv (Dun¬ 
lop et al. 2011, 2013; Wilkins et al. 2011; Finkelstein et al. 
2012a; Bouwens et al. 2014b; Duncan et al. 2014; Rogers 
et al. 2014). Disagreement between past observations on the 
existence or steepness of a color-magnitude relation (cf. Dun¬ 
lop et al. (2011) and Bouwens et al. (2011)) have recently 
been reconciled by Bouwens et al. (2014b) after addressing 
systematics in the selection and photometry between dif¬ 
ferent studies. Bouwens et al. (2014b) find a clear colour- 
magnitude relation (CMR) with bluer UV-slopes at lower 
luminosities, the relation is also found to evolve with bluer 
/3’s at high redshift (blue circles in Fig. 2). The existence of 
a strong colour-magnitude relation has also been confirmed 
by Rogers et al. (2014) at 2 ~ 5 for a sample of even greater 
dynamic range (pink diamonds in Fig. 2), measuring a CMR 
slope and intercept within error of the measured z ~ 5 val¬ 
ues of Bouwens et al. (2014b). 

While the evidence for a colour-magnitude relation at 
« ^ 6 is quite strong, during the crucial period of reion¬ 
ization at z ~ 7 the limited dynamic range and samples 
sizes of existing observations means that similar conclu¬ 
sions are less obvious. Given the importance of this period 
in reionization and also the importance the assumption of 
a colour-magnitude relation has on the conclusions of this 
study, it is pertinent to critically assess which model is sta¬ 
tistically favoured by the existing observations. To assess 
the statistical evidence for a colour-magnitude relation at 
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Figure 2. Observed average values of the UV continuum slopes 
/3 as a function of rest-frame UV magnitude, Muv, from Wilkins 
et al. (2011), Dunlop et al. (2011, 2013), Finkelstein et al. (2012a), 
Bouwens et al. (2014b), Duncan et al. (2014) and Rogers et al. 
(2014) at redshifts z ~ 4, 5, 6, 7 and 8 — 9. In the bottom panel, 
filled symbols show the average for 8 ~ 7 samples while the open 
symbol shows the averages for 2^9 (see respective papers for 
sample details and redshift ranges). 


2 ~ 7, we calculate the Bayesian information criterion (BIC, 
Schwarz 1978), defined as 

BIC = — 2 In £ max + k In N, (4) 

where £ max is the maximised value of the likelihood func¬ 
tion for the model in question (—2 In £ max = Xmin under the 
assumption of gaussian errors), k the number of parameters 
in said model and N the number of data-points being fit¬ 
ted. Values of ABIC greater than 2 are positive evidence 
against the model with higher BIC, whilst values greater 
than 6 (10) are strong (very strong) evidence against. The 
model fits were performed using the MCMC implementation 
of Foreman-Mackey et al. (2013) assuming a flat-prior. The 
resulting BIC for the constant /3 and linear Muv-dependent 
models are shown in Table 1. For the z ~ 7 observations, a 
color-magnitude relation is strongly favoured over a constant 

3 with a best-fit model of —2.05 ±0.04 — 0.13 ± 0.04(Muv + 

19.5). At z ~ 8, while the assumption of a constant 3 pro- 


Table 1. Bayesian Information Criterion (BIC) for the assump¬ 
tion of either a colour-magnitude relation or a constant /3, ABIC 
is defined as BlCconst —BIC s i ope . For each dataset, we also show 
the best-fit model parameters and corresponding 1 -a errors for 
the model with lowest BIC. 

Redshift BIC s i ope BlCconst ABIC 

2 ~ 7 35.2 a 45.0 9.8 

2~ 8 5.8 4.2 b -1.6 


“ /3(M UV ) = -2.05 ± 0.04 - 0.13 ± 0.04 x (M uv + 19.5) 
b /3= —2.00 ± 0.11 


vides a better fit (/3 = —2.00 ± 0.11), no model is strongly 
preferred over the other. 

Past studies into the ionizing emissivity during the EoR 
have often used a fixed average 3 = —2 to motivate or 
constrain £i on , e.g. Bolton & Haehnelt (2007), Ouchi et al. 
(2009), Robertson et al. (2010, 2013) and Kuhlen & Faucher- 
Giguere (2012). Although we now find good evidence for a 
colour-magnitude relation during this epoch, it does not nec¬ 
essarily make the assumption of a constant 3 invalid, as we 
must take into account the colours of the galaxies which 
dominate the SFR or luminosity density. In order to esti¬ 
mate an average 3 which takes into account the correspond¬ 
ing number densities and galaxy luminosities, we calculate 
(j9) , the average 3 weighted by the contribution to the 

total UV luminosity density: 


J7° ±uv(m) x <j>{m ) x (3) (m) 
IP Luv(rn) x Mm) 

° ^min 


(5) 


where Luv(ur), 4>{m) and (/3) (m) are the luminosity, num¬ 
ber density and average j3 at the rest-frame UV magnitude, 
m , respectively. We choose a lower limit of Luv = Muv = 
— 17, corresponding to the approximate limiting magnitude 
of the deepest surveys at z ^ 6. For the discrete bins in 
which (3) (m) is calculated, (/3) , this becomes a sum over 
the bins of absolute magnitude, k, brighter than our lower 
limit Muv = —17: 


< 0 ™ 


Y.Ik V,k X 4>k X (P) k 

Yk LuV,k X (j> k . 


( 6 ) 


The number density for a given rest-frame magnitude bin, 
( t>k , is calculated from the best-fitting UV luminosity func¬ 
tions of Bouwens et al. (2014a) at the corresponding redshift. 
We use the same luminosity function at each redshift for all 
of the observations for consistency. We note that given the 
relatively good agreement between estimates given their er¬ 
rors, the use of differing luminosity function estimates would 
have minimal effect on the calculated values. 

We estimate errors on {/3) uv through a simple 
Monte Carlo simulation, whereby (/3) m and the best-fitting 
Schechter (1976) parameters used to calculate 4>m are per¬ 
turbed by the quoted errors (making use of the full covari¬ 
ance measured by Bouwens et al. 2014a), this is repeated 
10 2 3 4 times. (3) and error are then taken as the median 
and lcr range of the resulting distribution. 
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Table 2. Bayesian Information Criterion (BIC) for the assump¬ 
tion of a redshift dependent or constant (/3)^ uv , ABIC is defined 
as BICconst— BIC 2 . For each dataset, we also show the best-fit 
model parameters and corresponding 1-cr errors for the model 
with lowest BIC. 

Redshift range BIC Z BICconst ABIC 

4 < ^ < 8 44.0 a 76.9 32.8 

5 < 2 < 8 41.0 b 46.2 5.2 

a (P) Pvv (z) = -1-59 ± 0.05 - 0.07 ± 0.01 x 2 . 
b Wpuv ( 2 ) = _L61 ± °- 12 “ °- 07 ± 0.02 x 2 - 
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Figure 3. Luminosity-weighted average p, (P) puv , as a function 
of redshift for the Muv — P observations shown in Fig. 2. The grey 
shaded region covers the range —2.2 < ft < —1.8 and the blue 
long-dashed line shows our parametrisation of (P) pvv vs z based 
on the observations of Bouwens et al. (2014b) (see Equation 7). 

Fig. 3 shows the calculated (P) Pvv and corresponding 
errors for each of the samples shown in Fig. 2. Overall, we 
can see that /3 ~ —2 is still a valid choice for a fiducial value 
of /3 during the epoch of reionization (z > 6) based on all 
of the existing observations. However, for several of the sets 
of observations there is evidence for a potential evolution at 
z < 7 (Wilkins et al. 2011; Finkelstein et al. 2012b; Bouwens 
et al. 2014b; Duncan et al. 2014). 

To determine whether there is statistical evidence for a 
redshift evolution in {P) Pvv , we again calculate the corre¬ 
sponding Bayesian information criteria to assess the relative 
merits of each model, assuming a simple linear evolution 
for the redshift dependent model. The resulting BIC and 
ABIC are listed in Table 2. Based on the observations in 
the redshift range 4 < 2 < 8, there is very strong evidence 
for a redshift-dependent (0) p over one which is constant. 
Because the model fits may be dominated by the 2 ~ 4 
observations and their significantly smaller errors, we also 
calculate the fits using only the 5 < 2 < 8 data. At z > 5, 
the collective observations do still favour redshift evolution, 
but the statistical significance is notably reduced. 

Throughout this work we choose to use the P observa¬ 
tions of Bouwens et al. (2014b) for both p vs Muv and (/3} uv 
as the basis of our analysis. This choice is motivated by the 
fact that these observations represent the largest samples 
studied (both in number and dynamic range) and due to 


the careful minimisation of the potential systematic errors 
are likely the least-biased observations. For this set of ob¬ 
servations, there is strong evidence for evolution in (/3) uv 
at 2 < 7. Specifically, from 2 ~ 4 to 2 ~ 7, (/3} uv steepens 
considerably from —1.9 ± 0.02 to —2.21 ± 0.14. Parametris¬ 
ing the Bouwens et al. (2014b) observations with a simple 
linear relation , we find: 

(P) pvv (2) = -1.52 ± 0.11 - 0.09 ± 0.02 x 2. (7) 

Based on this fit we predict a (P) Pvv = —2.24 for 2 ~ 8. 
Whilst this is significantly bluer than that based on the ex¬ 
isting observations at 2 ~ 8, it is comparable to the average 
P measured for fainter galaxies in the lower redshift samples 
and represents a reasonable extrapolation. We note that at 
2 s» 8, the (P) Puv changes significantly depending on the 
choice of average due to the significantly smaller samples 
observable and large scatter in the faintest bin. For exam¬ 
ple, using the bi-weight means recommended by Bouwens 
et al. (2014b), (P) Pvv = —1-74 based on their observations. 
Re-calculating using the inverse-weighted means of this same 
sample gives {P) Pl]v = —2.1, in better agreement with the 
observed trend at 2 < 8. 

While there is now good agreement on the existence and 
slope of the colour-magnitude relation between independent 
studies (cf. Bouwens et al. 2014b and Rogers et al. 2014), 
what is less well understood is the intrinsic scatter in the 
CMR and whether it is luminosity dependent. Currently, 
the most extensive study of the intrinsic scatter is that of 
Rogers et al. (2014), who found that the intrinsic scatter 
in P is significantly larger for bright galaxies. They also 
find an apparent lower limit (25th percentile) of P = —2.1 
which varies little with galaxy luminosity whilst the corre¬ 
sponding 75th percentiles increase significantly from fainter 
to brighter galaxies. Such a scenario implies that galaxies 
with P ^ —2.5 should be extremely rare at high-redshift, 
even though such galaxies are observed locally and at in¬ 
termediate redshifts (Stark et al. 2014b). Without a better 
understanding of the causes of the intrinsic scatter in P and 
the underlying stellar populations it is difficult to predict 
the expected numbers of such galaxies during this epoch. 

Rogers et al. interpret the intrinsic scatter as consistent 
with two simple scenarios: 1) the scatter is due to galaxy ori¬ 
entation, or 2) that brighter galaxies have more stochastic 
star-formation histories and the p variation is a result of ob¬ 
serving galaxies at different points in the duty cycle of star- 
formation. However, this second scenario is contrary to the 
theoretical predictions of Dayal et al. (2013), whereby fainter 
low-mass galaxies have more stochastic star-formation his¬ 
tories due to the greater effect of feedback shutting down 
star-formation in low-mass haloes. 

After the discovery that high-redshift galaxies exhibit 
significant UV emission lines by Stark et al. (2014a) (specif¬ 
ically Cm] at 1909A), it is worth asking if the presence of 
such far-UV emission lines can systematically affect mea¬ 
surements of the UV slope to the same degree which optical 
emission lines can affect age estimates and stellar masses. In 
Stark et al. (2014b), the authors find that the fitting of p is 
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not affected by the presence of UV emission lines in a sam¬ 
ple of young low-mass galaxies at 2 ~ 2. The same is true 
for galaxies out to z < 6 , where P is typically measured by 
fitting a power-law to three or more filters (Bouwens et al. 
2014b). However, at z ^ 7 where p must be measured using 
a single colour, the effect of UV emission line contamination 
in one of the filters is more significant. For a Cm] equiva¬ 
lent width of 13.5A (the highest observed in the Stark et al. 
(2014b) sample at z ~ 2) could result in a measured P which 
is too red by A/3 ~ 0.18 relative to the intrinsic slope based 
on the method outlined in Section 3. 

Given the limited samples of z > 6 galaxies with UV 
emission line detections, fully quantifying the effects of the 
emission line contamination on /3 observations is not possible 
at this time. As such, in this work we do not include UV 
emission lines in our simulated SEDs or attempt to correct 
for their effects on the observed /3s in Fig. 3. We do caution 
that despite the significant improvement on /3 measurements 
at high redshift, there may still be unquantified systematics 
when interpreting the UV slope during the EoR. 


3 MODELLING /3, £, ion AND n ion 

To model the apparent (3 ’s and corresponding emissivity 
coefficients for our two Lyman continuum escape models, 
we make use of composite stellar population models from 
Bruzual & Chariot (2003) (BC03). Using a stellar popula¬ 
tion synthesis code, we are able to calculate the full spectral 
energy distribution for a stellar population of the desired 
age, star-formation history, metallicity and dust extinction. 
Our code allows for any single-parameter star-formation his¬ 
tory (e.g. exponential decline, power-law, truncated or ‘de¬ 
layed’ star-formation models), and a range of dust extinction 
models. The models also allow for the inclusion of nebular 
emission (both line and continuum emission) proportional 
to the LyC photon rate, full details of which can be found 
in Duncan et al. (2014). 

For each resulting SED with known star-formation rate 
(SFR), we calculate the UV luminosity by convolving the 
SED with a top-hat filter of width 100A centred around 
1500A, as is standard practice for such studies (e.g. Finkel- 
stein et al. 2012b, McLure et al. 2013). To measure f3, each 
SED is redshifted to 2 ~ 7 and convolved with the WFC3 
F125W and F160W filter responses (hereafter J 125 and Hieo 
respectively). We then calculate /3 as: 

(3 = 4.43(Ji25 — Hieo) — 2. (8) 

as in Dunlop et al. (2013). This method is directly compa¬ 
rable to how the majority of the high redshift observations 
were made and should allow for direct comparison when in¬ 
terpreting the observations with these models. Using com¬ 
parable colours at z ~ 5 and 2 ~ 6 or different combinations 
of filters has minimal systematic effect on the calculated val¬ 
ues of P (see Dunlop et al. 2011 and Appendix of Bouwens 
et al. 2012 b). 

The LyC flux from these models is calculated before and 
after the applied absorption by gas (for nebular emission) 


and dust. We are therefore able to calculate the total escape 
fraction of LyC photons for a given stellar population. Using 
these values, it is therefore relatively straight-forward to link 
the observed P distribution of high-redshift galaxies with the 
distribution of predicted Ki on or Luv per unit SFR and the 
corresponding f eac ,tot. 

For the ionization-bounded nebula with holes model 
(Model A), the ‘observed’ SED is a weighted (proportional 
to /esc) sum of the un-attenuated starlight escaping through 
holes and the attenuated starlight and nebular emission 
from the H 11 and dust enclosed region. The two SED com¬ 
ponents are weighted proportional to the covering fraction 
(= 1 — /esc) of the H 11 and dust region. The resulting observ¬ 
able quantities are effectively the average over all possible 
viewing angles, as would be expected for a large sample of 
randomly aligned galaxies. 

In the case of Model B, the density bounded nebula, Ly¬ 
man continuum emission from the underlying stellar spec¬ 
trum is partially absorbed by the surrounding truncated 
Stromgren sphere with an escape fraction f esc ,neb- The re¬ 
maining Lyman continuum photons along with the UV- 
optical starlight and nebular emission are then attenuated 
by the surrounding dust shell according to the chosen dust 
attenuation law. The total escape fraction for this model is 
therefore 

f _ iq0.4 xA(LyC) r (Q\ 

J esc — ,/esc,neb V.' 7 / 

where A(LyC) is the magnitude of dust extinction for Ly¬ 
man continuum photons and is highly dependent on the 
choice of attenuation curve (see Section 3.1.6). 

3.1 Modelling assumptions: current constraints 
on stellar populations at z > 3 

Although there are now good constraints on both the UV 
luminosity function and UV continuum slope at high red¬ 
shift, both of these values suffer strong degeneracies with 
respect to many stellar population parameters. As such, con¬ 
straints on /esc, Cion and Kion still requires some assumptions 
or plausible limits set on the range of some parameters. In 
this section we outline the existing constraints on the rel¬ 
evant stellar population properties at high-redshift, discuss 
what assumptions we make in our subsequent analysis, and 
explore the systematic effects of variations in these assump¬ 
tions. 

3.1.1 Star-formation history 

Typically, star formation histories are parameterised as ex¬ 
ponential models, 

Psfr(I) oc exp(--), (10) 

r 

where r is the characteristic timescale and can be negative 
or positive (for increasing or decreasing SFH respectively). 
Or, alternatively as a power-law, following 

Psfr(1) oc (t — t() a . (11) 
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Figure 4. Power-law (blue continuous: Salmon et al. 2014, red 
dotted: Papovich et al. 2011) and exponential (green dashed: Pa- 
povich et al. 2011) fits to the median observed SFR-densities at 
z > 4 for 3 different star formation histories. The grey datapoints 
are taken from the compilation of SFR-density observations in 
Behroozi et al. (2013) with additional points from recent obser¬ 
vations (Smit et al. 2012; Bouwens et al. 2014a; Duncan et al. 
2014). For the power-law fits, the shaded red and blue regions 
correspond to the 1 -a errors on the slope of the power-law, a, 
quoted in the respective papers. The grey dot-dashed line shows 
the best-fit to the median star-formation rates across the the full 
cosmic history for the functional form outlined in Behroozi et al. 
(2013). All of the models assume an initial onset of star-formation 
at Zf = 12. 

The star-formation history (SFH) of high-redshift galaxies, 
SFR(t), is still very poorly constrained due to the limited 
rest-frame wavelengths available for SED fitting or spec¬ 
troscopy. For large samples of both intermediate and high 
redshift galaxies, it has been found that rising SFHs pro¬ 
duce better SED fits to the observed photometry (Maraston 
et al. 2010; Lee et al. 2014). However, reliably constraining 
the characteristic timescales, r or a, for individual galaxies 
at z > 2 is not possible for all but the brightest sources. 

Using a comoving number density selected sample of 
galaxies at high redshift, Papovich et al. (2011) found the 
average star formation history between 3 < 2 < 8 to be best- 
fitted by a power-law with a = 1.7± 0.2 or an exponentially 
rising history with r = 420 Myr. Recently, for the deep 
observations of the CANDELS GOODS South field, Salmon 
et al. (2014) applied an improved version of this method 
(incorporating the predicted effects of merger rates on the 
comoving sample) and found a shallower power-law with a = 
1.4 ± 0.1 produced the closest match. 

In Fig. 4, we show that all three of these models (a = 
1.4/1.7 and r = —450Myr) can provide a good fit to the ob¬ 
served evolution in the cosmic star-formation rate density at 
2 > 3 (Age of the Universe < 2 Gyr) through a simple scal¬ 
ing alone. However, the power-law fit with a = 1.4 provides 
the best fit to not only to the evolution of the SFR-density 
at ages < 2 Gyr, but also to the SFR density at later epochs. 
A smoothly rising star-formation is also favoured by hydro¬ 


dynamic models such as Finlator et al. (2011) and Dayal 
et al. (2013), although the star-formation histories of in¬ 
dividual galaxies are likely to be more varied or stochastic 
(Dayal et al. 2013; Kimm & Gen 2014). Furthermore, there is 
also growing evidence of galaxy populations with older pop¬ 
ulations and possibly quiescent populations (Nayyeri et al. 
2014; Spitler et al. 2014) suggesting some galaxies form very 
rapidly at high-redshift before becoming quenched. The as¬ 
sumption of a single parametrised SFH is clearly not ideal, 
however our choice of a rising power-law SFH with a = 1.4 is 
at least well motivated by observations and a more physical 
choice than a constant or exponentially declining SFH. 

3.1.2 Initial mass function 

Interpretation of extragalactic observations through mod¬ 
elling and SED fitting is typically done assuming a universal 
bi-modal Milky Way-like initial mass function (IMF) such 
as Kroupa (2001)/Chabrier (2003) or the unimodal Salpeter 
(1955) IMF. However, there is now growing evidence for sys¬ 
tematic variation in the IMF of both nearby (van Dokkum 
& Conroy 2010; Treu et al. 2010; Cappellari et al. 2012; Con¬ 
roy & van Dokkum 2012; Ferreras et al. 2013) and distant 
(Martfn-Navarro et al. 2014) early-type galaxies. 

Under a hierarchical model of galaxy evolution with 
downsizing, the bright galaxies in overdense regions observed 
at 2 > 3 are likely to eventually form into the massive early- 
type galaxies in which these variations can be found. Varia¬ 
tion in the slope of the IMF would have a significant effect 
on many of the critical observable properties at high redshift 
such as stellar masses and mass to UV light ratio’s. How¬ 
ever, given the lack of theoretical understanding as to how 
the IMF should vary with physical conditions, incorporat¬ 
ing the effects of a varying IMF at high-redshift is beyond 
the scope of this work. Throughout the following analysis 
we use the Chabrier (2003) IMF as our primary assumption, 
but also consider the systematic effect of a steeper IMF such 
as Salpeter (1955) on the inferred values or observables in 
Appendix A. 

3.1.3 Metallicity 

Current spectroscopic constrains on galaxy metallicities at 
2^3 indicate moderately sub-solar stellar and gas-phase 
metallicities (Shapley et al. 2003; Maiolino et al. 2008; 
Laskar et al. 2011; Jones et al. 2012; Sommariva et al. 2012). 
In addition, Troncoso et al. (2014), present measurements 
for 40 galaxies at 3 < 2 < 5 for which the observed metallic¬ 
ities are consistent with a downward evolution in the mass- 
metallicity relation (with increasing redshift). 

Measurements of galaxy metallicities at higher redshift 
(2 > 6) are even fewer due to the lack of high-resolution 
rest-frame optical spectroscopy normally required to con¬ 
strain metallicity. However, thanks to a clear detection of 
the Cm] emission line (1909A) and strong photometric con¬ 
straints, Stark et al. (2014a) are able to measure a metal¬ 
licity of « l/20th solar metallicity for a lensed galaxy at 
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2 = 6.029. Given these observations and the metallicities 
available in the Bruzual & Chariot (2003) models, we as¬ 
sume a fiducial metallicity of sa l/5th solar metallicity 
(Z = 0.004 = 0.2Z o ). 

3.1.4 Age 

For the rising star-formation history used throughout this 
work, there is a weak evolution of /3 as a function of age 
(A/3 « 0.13 between t « 100 Myr and « 1 Gyr) whereby 
older stellar populations have redder UV continuum slopes. 
However, at very young ages, the contribution of nebular 
continuum emission in the UV continuum can also signifi¬ 
cantly redden the apparent /3 compared to the much steeper 
underlying intrinsic UV slope (Robertson et al. 2010). This 
results in a degeneracy with respect to /3 between young and 
old populations. For example, for identical observed (stellar 
+ nebular continuum) UV-slopes, £i on (and Ki on ) for a young 
stellar population can be a factor of up ~ 0.5 (0.2) dex 
higher. Given this degeneracy, additional constraints from 
other parts of the electromagnetic spectrum are required in 
order to make a well informed choice of stellar population 
age. 

Due to the observational restrictions on high-resolution 
rest-frame optical spectroscopy, measurement of stellar pop¬ 
ulation ages for high-z galaxies is limited to photometric fit¬ 
ting and colour analysis. At a ~ 4, where the Balmer break 
can be constrained through deep Spitzer IR AC photometry, 
estimates of the average stellar population ages vary sig¬ 
nificantly from ~ 200 — 400 Myr (Lee et al. 2011) to ~ 1 
Gyr (Oesch et al. 2013) (dependent on assumptions of star- 
formation history). 

For galaxies closer to the epoch of reionization, con¬ 
straints on the Balmer/D(4000) breaks become poorer due 
to the fewer bands available to measure the continuum above 
the break, a problem which is exacerbated by the additional 
degeneracy of strong nebular emission lines redshifted into 
those filters (Schaerer & de Barros 2009, 2010). The effect 
of incorporating the effects of emission lines on SED fits at 
z ^ 5 is that on average the best-fitting ages and stellar 
masses are lowered. This is because the rest-frame optical 
colours can often be well fit by either a strong Balmer break 
or by a significantly younger stellar population with very 
high equivalent width Ha or Om emission. 

Recent observations of galaxies at high-redshift with 
constraints on the UV emission-line strengths (Lya or oth¬ 
erwise) have found that single-component star-formation 
histories are unable to adequately fit both the strong line 
emission and the photometry at longer wavelengths (Ro¬ 
driguez Espinosa et al. 2014; Stark et al. 2014a). For ex¬ 
ample, in order to match both the observed photometry 
at rest-frame optical wavelengths and the high-equivalent 
width UV emission lines of a lensed galaxy at 0 = 6.02, 
Stark et al. (2014a) require two stellar populations. In com¬ 
bination with a ‘young’ 10 Myr old starburst component, 
the older stellar component is best fitted with an age since 
the onset of star-formation of sa 500 Myr. 


Based on the observations discussed above and the me¬ 
dian best-fit ages found by Curtis-Lake et al. (2013), the 
star-formation history outlined in Section 3.1.1 is consistent 
with the existing limited constraints. At z ~ 7, the redshift 
of strongest interest to current studies of reionization, the 
assumed onset of star-formation at 2 = 12 gives rise to a 
stellar population age (since the onset of star-formation) of 
~ 390 Myr. 

3.1.5 Nebular continuum and line emission 

If nebular line emission is “ubiquitous” at high-redshift as an 
increasing number of studies claim (e.g. Shim et al. (2011); 
Stark et al. (2013); Smit et al. (2013)), the accompanying 
nebular continuum emission should also have a strong effect 
on the observed SEDs of high-redshift galaxies (Reines et al. 
2009). In this work, we include both nebular continuum and 
optical line emission using the prescription outlined in Dun¬ 
can et al. (2014) (and equivalent to the methods described 
in Ono et al. (2010); Schaerer & de Barros (2010); McLure 
et al. ( 2011 )). 

Whilst the strength of nebular emission in this model 
is directly proportional to the number of ionizing photons 
produced by the underlying stellar population, additionally 
both the strength and spectral shape of the nebular contin¬ 
uum emission are also dependent on the continuum emission 
coemcient, 7 ^ , given by 

(total) (HI) (2g) ( ge /) n(He + ) (Hell) n(He + + ) 

1V IV I IV I ( TT4-\ ' l u t TT-i-\ 

n{H + ) n(H + ) 

( 12 ) 

where 7 i HI \ 7 i HeI \ yl HeII ' > and y' 2q ' ) are the continuum 
emission coefficients for free-free and free-bound emission 
by Hydrogen, neutral Helium, singly ionized Helium and 
two-photon emission for Hydrogen respectively (Krueger 
et al. 1995). As in Duncan et al. (2014), the assumed con¬ 
tinuum coefficients are taken from Osterbrock & Ferland 
(2006), assuming an electron temperature T = 10 4 K and 
electron density n e = 10 2 cm -3 and abundance ratios of 
y + = » ( (g E +) ) = O - 1 and V ++ = n< a(H + + + , ) = 0 (Krueger et al. 
1995; Ono et al. 2010). 

Although the exact ISM conditions and abundances of 
high-redshift Hu regions is not well known, singly and dou¬ 
bly ionised helium abundances for nearby low-metallicity 
galaxies have been found to be y + « 0.08 and y ++ 0.001 

(Dinerstein & Shields 1986; Izotov et al. 1994; Hagele et al. 
2006). We estimate that for the age, metallicity and dust val¬ 
ues chosen for our fiducial model (see Table 3), variations 
of A y + = 0.05 corresponds to A/3 = 0.004, while values of 
y ++ as large as 3% (Izotov et al. 2013) would redden the 
observed UV slope by A/3 = +0.003. In this case, because 
the nebular continuum emission is dominated by the stel¬ 
lar continuum at these wavelengths for our assumption, the 
effects of variation in the Hu region properties is negligible 
and our interpretation of the observed UV slopes should not 
be affected by our assumed nebular emission properties. 
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3.1.6 Dust Extinction 

In Bouwens et al. (2009), it is argued that the most likely 
physical explanation for variation in the observed /3s be¬ 
tween galaxies is through the variation in dust content. 
Large changes in metallicity and ages are required to pro¬ 
duce significant variation in /3 (see later discussion in Sec¬ 
tion 3.2), however as previously discussed in this section, 
such large variations are not observed in the galaxy popu¬ 
lation in either age or metallicity at z > 3 based on current 
observations. For the fiducial model in our subsequent anal¬ 
ysis, we allow the magnitude of dust extinction (Ay) to vary 
along with / eBC , but we must also choose a dust attenuation 
law to apply. 

Direct measurements of dust and gas at extreme red- 
shifts are now possible thanks to the sub-mm facilities of 
ALMA. However, the current number of high redshift ob¬ 
servations is still very small. Ouchi et al. (2013) and Ota 
et al. (2014) observe Lyman alpha emitters (LAEs) at z ~ 7, 
finding only modest dust extinction ( E(B — V)~ 0.15). Re¬ 
cent work by Schaerer et al. (2014) extends the analysis to a 
larger sample of five galaxies, finding a range in dust extinc¬ 
tion of 0.1 < Ay < 0.8. While the dust attenuation in these 
objects is consistent with normal extragalactic attenuation 
curves such as (Calzetti et al. 2000) or the SMC extinction 
curve (e.g. Pei (1992)), the results are not strong enough to 
constrain or distinguish between these models. Similarly, for 
broadband SED fits of high-redshift objects, neither a star- 
burst or SMC-like attenuation curve is strongly favoured 
(Salmon et al. 2014). Based on these factors, we assume the 
starburst dust attenuation curve of Calzetti et al. (2000) in 
order to make consistent comparisons with the SED fitting 
of Duncan et al. (2014) and Meurer et al. (1999) dust correc¬ 
tions to UV star-formation rates (e.g. Bouwens et al. (2011); 
Smit et al. (2012)). 

In addition, due to the lack of constraints on the dust 
attenuation strengths at wavelengths less than 1200A we 
must also assume a plausible extrapolation below these 
wavelengths. For our fiducial model, we simply extrapo¬ 
late linearly based on the slope of the attenuation curve 
at 1200 - 1250A, in line with similar works on the escape 
fraction of galaxies (Siana et al. 2007). In addition, we also 
assume a second model in which the extreme-UV and Ly¬ 
man continuum extinction follows the functional form of 
the component of the Pei (1992) SMC extinction model at 
A ^ 1000A, whereby the relative absorption begins to de¬ 
crease below 800A. The systematic effect of choosing this 
second assumption along with a third assumption of an SMC 
extinction curve are outlined in Table A. 

As shown in Fig. 1, we assume a simple foreground dust 
screen (Calzetti et al. 1994) and that dust destruction is min¬ 
imal and/or balanced by grain production (Zafar & Watson 
2013; Rowlands et al. 2014); effectively that dust for a given 
model is fixed with relation to the stellar population age. 
The assumption of a different dust geometry, such as one 
with clouds dispersed throughout the ISM, would require a 
greater amount of dust to achieve the same optical depth 
and could also have a significant effect on the extinction of 


nebular emission relative to that of the stellar continuum 
(Zackrisson et al. 2013). 


3.1.7 Differing SSP models: the effects of stellar rotation 
and binarity 

The choice of Bruzual & Chariot (2003) stellar population 
synthesis (SPS) models in this work was motivated by the 
more direct comparison which can be made between this 
analysis and the stellar mass, luminosity and colour mea¬ 
surements based on the same models, e.g. Finkelstein et al. 
(2012a); Duncan et al. (2014). However, several other SPS 
models are available and in common usage, e.g. Starburst99 
(Leitherer et al. 1999), Maraston (2005) and FSPS (Conroy 
et al. 2009a, b). 

Due to differences in assumptions/treatment of vari¬ 
ous ingredients such as horizontal branch morphology or 
thermally-pulsating asymptotic giant branch (TP-AGB) 
stars, the SEDs produced for the same input galaxy prop¬ 
erties (such as age and metallicity) can vary significantly. 
The full systematic effects of the different assumptions and 
models for galaxies at high-redshift is not well quantified and 
adequately doing so is beyond the scope of this work. We do 
however caution that these systematics could significantly 
affect the inferred ionizing photons rates of galaxies during 
the EoR. In particular, it has been suggested that rotation 
of massive stars could have a significant effect on the UV 
spectra and production rate of ionizing photons (Vazquez 
et al. 2007). 

In Leitherer et al. (2014), the effects of the new stellar 
models including stellar rotation (Ekstrom et al. 2012) are 
incorporated into the Starburst99 SPS models. The resulting 
SEDs are changed drastically with an increase in the ionizing 
photon rate of up to a factor of five for the most extreme 
model of rotation. 

A second, equally significant effect comes from the in¬ 
clusion of binary physics in stellar population synthesis mod¬ 
els. It is now believed that the majority of massive stars ex¬ 
ist in binaries (Sana et al. 2012, 2013; Aldoretta et al. 2014) 
while the majority of stellar population models (including 
all of the aforementioned SPS models) are for single stars. 
The BPASS code of Eldridge & Stanway (2009, 2011)) incor¬ 
porates the physics of binary rotation on massive stars to 
explore the effects on the predicted stellar population fea¬ 
tures, observing similar effects to the addition of rotation in 
single star models; an increased fraction of red supergiants 
which go on to form bluer UV bright Wolf-Rayet stars. 

We note that for the same assumed stellar populations 
parameters (age, metallicity, SFH, dust), the use of either 
of these models would result in SEDs with bluer UV contin¬ 
uum slopes and an increased LyC production rate. However, 
the full ramifications of how these models may change the 
interpretation of SEDs, stellar masses and /3s for the obser¬ 
vations of high-redshift galaxies is beyond the scope of this 
work. 
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Figure 5. Left: UV continuum slope [3 as a function of total escape fraction, /esc, and dust extinction, Ay, for the ionization bounded 
nebula with holes continuum escape model (Model A, Fig. 1 left) with stellar population properties as outlined in Section 3.1. The 
contours indicate lines of constant f3 around the observed average (3, and the light grey arrow indicates how those contours move for 
a stellar population with solar metallicity. Middle and right: log 10 fe sc£ion and log 10 f e sc^ion as a function of escape fraction and dust 
extinction respectively for the same continuum escape model. Solid contours represent lines of constant f e sc£ion and /esc^ion whilst the 
dashed contour shows where (3 = —2 is located for reference. The green labelled contour shows the assumed f esc £ion = 24.5 value of R13 
and the equivalent in /esc^ion ( see text for details). For the colour scales below the centre and right panels, the lower black tick labels 
correspond to the scale for the fiducial model ( Z = 0.2 Z@) whilst the grey upper tick label indicate how f e S c£ion and f e sc^ion change for 
stellar populations with solar metallicity. 
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Figure 6. Left: UV continuum slope (3 as a function of Hu region escape fraction, f es c , and dust extinction, Ay, for the density bounded 
nebula continuum escape model (Model B, Fig. 1 right) with stellar population properties as outlined in Section 3.1. The contours indicate 
lines of constant (3 around the observed average /3, and the light grey arrow indicates how those contours move for a stellar population 
with solar metallicity. Centre and right: log 10 fe sc£ion and log 10 /esc^ion (where f esc is the total dust attenuated escape fraction) as 
a function of escape fraction and dust extinction respectively for the same continuum escape model. Solid contours represent lines of 
constant / e sc Cion and /esc^ion whilst the dashed contour shows where (3 = — 2 is located for reference. The green labelled contour shows 
the assumed / e sc£ion — 24.5 value of R13 and the equivalent in /esc^ion (see text for details). For the colour scales below the centre 
and right panels, the lower black tick labels correspond to the scale for the fiducial model (Z = 0.2 Z@) whilst the grey upper tick label 
indicate how f esc £ion and /esc^ion change for stellar populations with solar metallicity. 
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3.2 Observed UV slopes as a function of / eac and 
dust extinction 


Table 3. Summary of the stellar population assumptions for our 
fiducial 0 = —2 model. 


The galaxy properties with the largest uncertainties and ex¬ 
pected variation are the optical depth of the dust attenua¬ 
tion (or extinction magnitude Ay) and the property we wish 
to constrain photometrically, the escape fraction of ionizing 
photons /esc- For the assumptions of our fiducial model (Ta¬ 
ble 3), we want to explore the possible range of these two 
properties which are consistent with the observed UV slopes 
(as calculated for each model following Eq. 8) and what con¬ 
straints can then be placed on the ionizing emissivity coef¬ 
ficients /esc/ion Or /esc^ion- 

In the left panels of Fig. 5 and Fig. 6, we show how /3 
varies as a function of the dust extinction magnitude and es¬ 
cape fraction for each of the continuum escape mechanisms 
respectively (Section 2/Fig. 1). For both models, 0 is rela¬ 
tively constant as a function of f esc at low values of escape 
fraction (/ esc < 0.1) . At larger escape fractions, the two 
mechanisms produce different UV slopes. For Model A, as 
/esc increases to ~ 10% (covering fraction « 90%) the unat¬ 
tenuated stellar continuum begins to dominate the overall 
colours as / eS c increases and by f esc > 30% the observed 
average 0 is determined only by the unattenuated light, ir¬ 
respective of the magnitude of the dust extinction in the 
covered fraction. This effect is also illustrated in a different 
way in Figure 9 of Zackrisson et al. (2013), whereby the same 
amount of dust extinction in the covered/high-density re¬ 
gions has a decreasing effect on the observed 0 as the escape 
fraction increases. This means that if Lyman continuum is 
escaping through holes in the ISM and is un-attenuated by 
dust, it is possible set constraints on the maximum escape 
fraction possible which is still consistent with the UV slopes 
observed. 

For Model B, 0 remains constant with f eS c,neb at a fixed 
dust extinction until f esc « 30%. Beyond this, the reduction 
in nebular continuum emission from the high escape fraction 
begins to make the observed 0s bluer for the same magnitude 
of dust extinction. 

For both escape mechanisms, a UV slope of 0 ~ —2 is 
achievable with only moderate amounts of dust extinction 
required {Ay « 0.4 and 0.25 for models A and B respectively 
at /esc ~ 0.2). When metallicity is increased to Z = Z Q , the 
isochromes (of constant 0) are shifted downwards such that 
0 ~ — 2 requires negligible dust attenuation (cf. Robertson 
et al. 2013). In Section 3.4 we further explore the effects 
of varying metallicity on the apparent 0 and corresponding 
emissivity coefficients. However, first we wish to examine the 
range of £i orl and Ki on which correspond to the values of f eBC 
and Ay consistent with 0 ~ — 2 found here. 


3.3 £i on and Ki on as a function of f eBC and dust 
extinction 

In the centre and right panels of Figures 5 and 6 we show 
logio /escCion and log 10 /escKion as a function of /esc and the 
extinction magnitude of dust in the covered regions {Ay). 
For dust model A, the ionization bounded nebula with 


Star-formation history 
Initial Mass Function 
Dust attenuation curve 
Metallicity 
Nebular Emission 
Age 


SFR oc t 1A a 
Chabrier (2003) 
Calzetti et al. (2000) 
Z = O.2Z 0 

Continuum included^ 
390 Myr c 



Model A 

Model B 

Dust attenuation magnitude Ay 

0.31 

0.27 

Escape fraction f eS c,neb 

0.16 

0.42 

logio fe sc£ion 

24. 5 d 

24.5 

log 10 /esc^ion 

52.39 

52.34 


a Salmon et al. (2014) 

k T = 10 4 K, n e = 10 2 cm -3 , y + = 0.1 and y ++ = 0 
c Age at z ~ 7 since onset of star-formation at z ~ 12. 
d The assumed log 10 / es c£ion of R13 


holes, there is very little dependence of the ionizing pho¬ 
ton rate per unit UV luminosity on the magnitude of dust 
extinction in the covered fraction (centre panel). Because 
the increase dust extinction magnitude around the high col¬ 
umn density areas only affects the UV/optical component, 
the increasing dust extinction results in higher values of 
logic /esc£ion due to the increased absorption of the UV light 
in the dust covered regions. For this dust model, the corre¬ 
sponding ionizing photon rate per unit SFR (/ es cftion) has 
zero evolution as a function of dust extinction in this geom¬ 
etry. 

The assumed log 10 / esc £ion = 24.5 of R13 is consistent 
with 0 = —2 for this model, with an escape fraction of 
/esc = 0.16 and moderate dust extinction {Ay = 0.31). 
Given the low escape fractions which are still consistent with 
blue 0 slopes, a value of log 10 / esc £ion = 24.5 does represent 
a relatively optimistic assumption on the ionizing efficiency 
of galaxies. However, it is still rts 0.2 dex lower than the 
largest f eBC still consistent with a UV slope of 0 = —2. 

In contrast to model A, because the dust in the density 
bounded nebula (model B) is assumed to cover all angles, 
increases in the dust extinction magnitude results in signifi¬ 
cantly smaller log 10 / eS c£ion/log 10 / esc K ion for the same fixed 
/esc- This can be seen clearly in the centre and right panels 
of Fig. 6. 

Despite this, an assumed value of log 10 f e sc £ion = 24.5 
is still consistent with 0 = —2 for this model. However, it re¬ 
quires a higher escape fraction (/ eB c = 0.42) and lower dust 
extinction {Ay = 0.27) to achieve this for the same under¬ 
lying stellar population. The maximum f esc ,neb = 1 is still 
consistent with the fiducial UV slope, but the increased dust 
required to match 0 = —2 means that the total LyC escape 
fraction is reduced and that the corresponding maximum 
log 10 /esc£ion is only marginally higher than the assumptions 
of (Kuhlen & Faucher-Giguere 2012) or (R13). 

For the assumed stellar population properties in our 
reference model, the UV continuum slope of the intrinsic 
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dust-free stellar population is /3 = —2.55 excluding the con¬ 
tribution of nebular continuum emission. This value is sig¬ 
nificantly bluer than the dust-free /3 assumed by the Meurer 
et al. (1999) relation commonly used to correct UV star- 
formation rates for dust absorption. It is however in better 
agreement with the dust-free /3’s estimated recently for ob¬ 
served galaxies at z ^ 3 (Castellano et al. 2014; de Barros 
et al. 2014) and the theoretical predictions of Dayal & Fer¬ 
rara (2012). 

3.4 Effect of different stellar population 
properties on £i on and Ki on vs /3 

Given the strong evidence for both luminosity and redshift 
dependent /3s, we wish to explore whether evolution in each 
of the stellar population parameters can account for the ob¬ 
served range of /3s and estimate what effect such evolution 
would have on the inferred values of / es c£ion and / e scKi on . 
At 2 > 6, the the average /3 for the brightest and faintest 
galaxies by A/3 « 0.6 (Fig. 2). As a function of redshift, 
the evolution in /3 is less dramatic, with average slopes (at 
a fixed luminosity) reddening by A/3 ~ 0.1 in the ~ 400 
million years between z ~ 7 and z ~ 5. 

In the left panel of Fig. 7, we show how /3 and / eS c£ion 
or /esc Ki on vary as a function of each model parameter for 
the ionization-bounded nebula model (Model A) with the re¬ 
maining parameters kept fixed at our fiducial /3 = —2 model 
(Table 3). The corresponding values for the density bounded 
nebula model (Model B) are shown in the right-hand panel 
of Fig. 7. How / esc £ion and / es cKi on vary with respect to /3 
for changes in the different parameters differs significantly: 

• Dust: As has already been seen in Fig. 5, /3 varies 
strongly as a function of dust attenuation strength for the 
ionization-bounded nebula model. There is however minimal 
variation in / esc £ion or / es cKion with respect to that large 
change in /3. The inferred / esc £ion or / eS cKi on can justifiably 
be considered constant as a function of redshift for this es¬ 
cape model if it is assumed that evolution in the dust ex¬ 
tinction is responsible for the observed evolution of /3. 

For model B, the density-bounded nebula, the large evo¬ 
lution in /3 is coupled to a significant evolution in the in¬ 
ferred emissivity coefficients. For a change in the UV slope 
of A/3 « 0.2, there is a corresponding evolution in / eS c£ion 
and/escKkm of « 0.19 and 0.25 dex respectively. 

• Metallicity: Between extremely sub-solar (Z = 
0.005 Zq) and super-solar (Z > 1 Zq) metallicities, the 
variation in /3 is large enough to account for the wide range 
of observed average /3’s for both Lyman escape mechanisms. 
In this regard, metallicity evolution is a plausible mecha¬ 
nism to explain the apparent variation in /3. However, such 
a wide variation in metallicities is not supported by current 
observations (see Section 3.1.3) . 

Both / e action and /escKion vary by a factor of ~ 2 across 
the full metallicity range modelled in this work, with bluer 
low-metallicity stellar populations producing more Lyman 
continuum photons per unit SFR/UV luminosity. 

• Age: Due to the strong nebular continuum contribution 


to the overall spectra at young ages, evolution in the stellar 
population age results in a more complicated /3-emissivity 
relation. When nebular emission is included the continuum 
emission reddens the slope at very young ages before turning 
over at t ~ 100 Myr and reddening with age towards ages 
of t = 1 Gyr and greater. For both LyC escape models, 
younger stellar populations results in a higher number of 
ionizing photons per unit UV luminosity. 

For Model A, the effect of reddening by nebular emission 
at young ages is more pronounced due to the lower nebular 
region escape fraction ( fesc,neb = 0.21) in our fiducial model. 
The result of this reddening is that for the same /3 = —2, the 
corresponding log 10 / esc £ion can be either « 24.5 or « 25.. 
This represents a potentially huge degeneracy if the ages of 
galaxy stellar populations are not well constrained. 

Between z ~ 7 and z ~ 4, the evolution in stellar popula¬ 
tion age can only account for a reddening of A/3 « 0.1, less 
than the evolution in both the normalisation of the colour- 
magnitude relation and in {P) pvv ■ Similarly, it worth noting 
that at z ~ 7, the assumption of a later onset for star- 
formation (e.g. z ~ 9, Planck Collaboration et al. (2015)) 
results in a ~ 0.14 dex increase in the intrinsic £i on and a 
change in the UV slope of A/3 ~ —0.065. 

• fe sc- For both models of LyC escape, variation in f esc 
has a strong evolution in f e sc £i on orf eBC Ki on with respect 
to changes in /3. However, for Model B (density bounded 
nebula) the range of /3 covered by the range of fesc,neb 
(0 ^ fesc,neb ^ 1) is only ~ 0.2 dex, significantly less than 
the range of colours reached by variation in the other stellar 
population parameters. 


4 ESTIMATED GALAXY EMISSIVITY 
DURING REIONIZATION 

Using our improved understanding of the ionizing efficiencies 
of galaxies during the EoR, we can now estimate the total 
ionizing emissivity N ; on of the galaxy population at high 
redshift following the prescription outlined in Equations 2 
and 3. We quantify puv and psfr using the latest available 
observations of the galaxy population extending deep in the 
epoch of reionization. 

4.1 Observations 

Thanks to the deep and wide near-infrared observations of 
the CANDELS survey (Grogin et al. 2011; Koekemoer et al. 
2011) and the extremely deep but narrow UDF12 survey 
(Koekemoer et al. 2013), there now exist direct constraints 
on the observed luminosity function deep into the epoch of 
reionization. In this paper, we will make use of the UV lu¬ 
minosity functions calculated by McLure et al. (2013) and 
Schenker et al. (2013) at z = 7 —9 as part of the UDF12 sur¬ 
vey along with the recent results of Bouwens et al. (2014a) 
and (Finkelstein et al. 2014) at z ^ 4 and the latest results 
from lensing clusters at z > 8 (Oesch et al. 2014; McLeod 
et al. 2014). 
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Model A: Model B: 

Ionization bounded nebula with holes Density bounded nebula 






Figure 7. Evolution of log 10 / es c£ion (top panels) and log 10 ,/esr^ion (bottom panels) vs /3 as a function of changes in the stellar 
population age (green dashed), metallicity (red dotted), dust extinction (yellow dot-dashed) and escape fraction (blue continuous) with 
the remaining parameters fixed to the fiducial values listed in Table 3. Values are plotted for both the ionization bounded nebula with 
holes (Model A; left panels) and density bounded nebula (Model B; right panels). Some individual points are labelled for both Lyman 
continuum escape models to illustrate the range and differences in evolution between each model. Note that the small difference in /3 
between Model A and Model B for the case of zero-dust (0 Ay) is due to the difference in nebular emission contribution for the two 
models; f esc ,neb = 0.16 and / eS c,ne6 = 0.42 for models A and B respectively. 


A second, complimentary constraint on the amount of 
star-formation at high redshift is the stellar mass function 
and the integrated stellar mass density observed in subse¬ 
quent epochs. As the time integral of all past star-formation, 
the stellar mass density can in principal be used to constrain 
the past star-formation rate if the star-formation history is 
known (Stark et al. 2007). A potential advantage of using 
the star-formation rate density in this manner is that by be¬ 
ing able to probe further down the mass function, it may be 
possible to indirectly measure more star-formation than is 
directly observable at higher redshifts. Or to outline in other 
terms, if the total stellar mass density of all galaxies can be 
well known at 2 ~ 4 or z ~ 5, strict upper limits can be 
placed on the amount of unobserved (e.g. below the limiting 
depth of 2 ~ 8 observations) or obscured star-formation at 
z > 6. 


At its simplest, the relation between a star-formation 
history, S(t), and the resulting stellar mass M* (or stellar 


mass density p*) is given by 

M,{t z ) = (1 - e z ) X [ ‘ S{t) dt (13) 

Jt f 

where t/ and t z are the age of the Universe at the onset of 
star-formation and observed redshift respectively and t z is 
the fraction of mass returned to the ISM at the observed red¬ 
shift. For a parametrised star-formation history, F(t), which 
is normalised such that F(t) dt = lM 0 (Mpc -3 ), we can 

substitute S(t) = C 0 bs,z x F(t). The normalisation, C 0 bs,z, 
accounts for the normalisation of the star-formation history 
required to match the observed stellar mass (or stellar mass 
density) at the redshift 2 . 

From recent observations of the stellar-mass density at 
2 > 4, such as those from the stellar mass functions pre¬ 
sented in Duncan et al. (2014) and Grazian et al. (2014), it is 
then straight-forward to calculate the corresponding C 0 b s ,z 
and the inferred star-formation history ( C 0 bs,z x F(t)) for 
any assumed parametrisation. Motivated by the discussion 
of star-formation histories in Section 3.1, we assume a nor- 
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malised star-formation history which is F(t) oc t 1 ' 4 (Salmon 
et al. 2014), as is done for the modelled (5 values. 

Using our stellar population models, we calculate e for 
the star-formation history and metallicity used in this work 
at any desired redshift. We find that for a Chabrier (2003) 
IMF, t z varies from e z « 0.29 at 2 = 8 to t z « 0.35 at 2 = 4 
(for a Salpeter IMF and constant SFH at ~ 10 Gyr old, we 
calculate e = 0.279 in agreement with the figure stated in 
R13). 

4.2 Mon for constant / esc £ ion and f esc K ion 

In Fig. 8, we show the estimated ionizing emissivity as a 
function of redshift for UV luminosity function observations 
of Bouwens et al. (2014a); Finkelstein et al. (2014); Oesch 
et al. (2014) and McLeod et al. (2014), assuming our fidu¬ 
cial constant log 10 / esc £ion = 24.5 (Table 3, as assumed in 
Robertson et al., 2013). The integrated luminosity density 
and corresponding confidence intervals for each UV LF ob¬ 
servation are estimated by drawing a set of LF parameters 
from the corresponding MCMC chain or likelihood distribu¬ 
tion obtained in their fitting. This is repeated 10 4 times to 
give a distribution from which we plot the median and 68% 
confidence interval. The luminosity functions of Bouwens 
et al. (2014a) and Finkelstein et al. (2014) span a large red¬ 
shift range from 2 ~ 4 to z ~ 8 and predominantly make use 
of the same imaging data (including the ultra deep UDF12 
observations, Koekemoer et al. 2013) but use different re¬ 
ductions of said data and differing selection and detection 
criteria, we refer the reader to the respective papers for more 
details. 

Also shown are the current constraints at 2 ~ 10 based 
on the luminosity function from Oesch et al. (2014). Due 
to the small number of sources available at 2 ~ 10 and 
the very large uncertainty in their redshift, the luminosity 
function is not well constrained at these redshifts. Fits us¬ 
ing the Schechter (1976) parametrisation realistically allows 
only one free parameter to be varied while the remainder 
are fixed to their 2 ~ 8 values (we plot the W on predicted 
for the luminosity function parameters where (/), is allowed 
to vary). Despite the large uncertainty, we include these val¬ 
ues to illustrate the early suggestions of Oesch et al. (2014) 
and other works (Zheng et al. 2012; Coe et al. 2012) that 
the luminosity function (and hence the inferred underlying 
star-formation rate density) begins to fall more rapidly at 
2^9 than extrapolations from lower redshift naively sug¬ 
gest. However, more recent analysis by McLeod et al. (2014) 
of existing Frontier Fields data is in better agreement with 
the predicted luminosity density at 2 « 9. Completed ob¬ 
servations of all six Frontier Fields clusters should provide 
significantly improved constraints at 2 > 9 (Coe et al. 2014), 
although the large cosmic variance of samples due to the 
volume effects of strong lensing will limit the constraints 
that can be placed at the highest redshifts (Robertson et al. 
2014). 

For all redshift samples plotted in Fig. 8, filled sym¬ 
bols correspond to the UV luminosity density integrated 


to Muv = —17 and Muv = —13 respectively. The pre¬ 
dicted Ni on for the UV luminosity functions of Schenker 
et al. (2013) and McLure et al. (2013) (not shown in Fig. 8) 
effectively reproduce the UV luminosity density constraints 
outlined in R13 and lie between those predicted by Bouwens 
et al. (2014a) and Finkelstein et al. (2014). The more re¬ 
cent works of Bouwens et al. (2014a) and Finkelstein et al. 
(2014) show a greater disagreement between both them¬ 
selves and previous works. This is a concern as it means 
that the choice of luminosity function (and hence the under¬ 
lying selection/methodology) could have a significant effect 
on the conclusions drawn on galaxies’ ability to complete or 
maintain reionization by the desired redshift. 

At 2 ~ 6, the UV luminosity density from galaxies 
brighter than the limiting depth observed by Bouwens et al. 
(2014a) is large enough to maintain reionization for a clump¬ 
ing factor of three. This is in contrast to the previous LF of 
Bouwens et al. (2012a) and the results of Finkelstein et al. 
(2014) which require a contribution from galaxies fainter 
than Muv = —17 (approximately the observational lim¬ 
its) to produce the Ni on needed to maintain reionization. 
At 2 ~ 8, the large uncertainties (and fitting degeneracies) 
in both the faint-end slope and characteristic luminosity of 
the luminosity function means that both of the luminosity 
density (and hence Ni on ) assuming a constant / esc £ ion) es¬ 
timates included in this work agree within their lcr errors. 
For the redshift dependent ionizing emissivity inferred by 
Bouwens et al. (2015 , assuming Chu = 3), faint galaxies 
down to at least Muv = —13 are required for all redshifts 
at 2 ^ 6 based on the current UV luminosity density con¬ 
straints and the fiducial / es c£km = 24.5. 

To convert the star-formation rates inferred by the stel¬ 
lar mass densities observed by Duncan et al. (2014) and 
Grazian et al. (2014) to an N[ on directly comparable with 
the LF estimates, we choose the f esc Kion at the log\of esc —Av 
values where /3 = —2 and log 10 / eS c£ion = 24.5. For the ion¬ 
ization bounded nebula with holes model, the corresponding 
logic /escKion = 52.39, whilst for the density bounded nebula 
logio fe sc^ion — 52.34. 

The ionizing photon rate predicted by the 2 ~ 4 stellar 
mass functions of Duncan et al. (2014) and Grazian et al. 
(2014) for stellar masses greater than 1O S ' 55 M0 (the es¬ 
timated lower limit from (Duncan et al. 2014) for which 
stellar masses can be reliably measured at 2 ~ 4) are 
shown as the blue and pink lines plotted in Fig. 8. Plot¬ 
ted as thick solid and dashed lines are the Ni on assuming 
logic /escKicm = 52.39 with the corresponding shaded area 
illustrating the systematic offset for logi 0 f e scKion = 52.34. 

The UV luminosity density at 2 > 4 implied by the 
SMD observations of Grazian et al. (2014) (for M > 
10 8 ' 55 Mq) are in good agreement with the UV LF esti¬ 
mates of Finkelstein et al. (2014) when integrated down to 
Muv = —17 (~ observation limits). However, when inte¬ 
grating the stellar mass function down to significantly lower 
masses such as > 10 7 Mq, the shallower low-mass slope of 
(Grazian et al. 2014) results in a negligible increase in the 
total stellar mass density (« 0.07 dex). This is potentially 
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Figure 8. Ionizing emissivity Ni on predicted for a fixed f e action for the UV luminosity function observations of Bouwens et al. (2014a), 
Oesch et al. (2014), Finkelstein et al. (2014) and McLeod et al. (2014). Filled symbols are for luminosity functions integrated down 
to Muv = —17 (typical limiting depth of while open symbols correspond to a UV luminosity density integrated down to a constant 
M\ jy = —13 for all redshifts. Also shown are the ionizing emissivities inferred by the z ~ 4 stellar mass density observations of Duncan 
et al. (2014) (solid blue line) and Grazian et al. (2014) (dashed pink line) for the power-law star-formation history outlined in Section 3.1. 
For the stellar mass function predictions, the thick solid and dashed line correspond to a constant log 10 /escKion = 52.39 (Lyman escape 
model A) while the filled region shows the systematic offset when assuming log 10 f e scKion = 52.34 for Lyman escape model B, see the 
text for details. We show the IGM emissivity measurements and corresponding total errors of Becker & Bolton (2013) (tan line and 
cross-hatched region respectively) and the ionizing background constraints of Bouwens et al. (2015) (blue-gray diagonal-hatched region). 
The ionizing emissivity required to maintain reionization as a function of redshift and clumping factors (Chu) of one, three and ten are 
shown by the wide grey regions (Madau et al. 1999, also see Equation 18 of Bolton & Haehnelt 2007). 


inconsistent with the star-formation (and resulting stellar 
mass density) implied at z > 6 when galaxies from below 
the current observations limits of the luminosity functions 
are taken into account (open plotted symbols). 

The higher stellar mass density observed by Duncan 
et al. (2014) results in Ni on most consistent with those in¬ 
ferred by the Bouwens et al. (2014a) luminosity density 
(Muv > —17). Overall, there is reasonable agreement be¬ 
tween the implied star-formation rates of the luminosity and 
stellar mass functions. However, both estimates have com¬ 
parable systematic differences between different sets of ob¬ 
servations. 

Since the observational limit for the stellar mass func¬ 
tion is currently limited by the reliability of stellar mass esti¬ 
mates for galaxies rather than their detection, better stellar 
mass estimates alone (through either deeper IRAC data or 
more informative fitting priors) could extend the observa¬ 
tional limit for current high redshift galaxy samples. Im¬ 
proved constraints on the stellar mass functions at z ^ 6 are 
therefore a viable way of improving the constraints on SFR 


density and ionizing emissivity of faint galaxies at higher 
redshifts. 


4.3 N ion for evolving / esc £km 

To estimate what effect a /? dependent / esc £ion or / es c/tion 
would have on the predicted M on , we assume two sepa¬ 
rate /esc£ion(,3) relations based on the predicted evolution 
of f e s C £ion and / esc «ion vs (3 shown in Fig. 7. The first rela¬ 
tion, ‘ModelB_dust’, is based on the relatively shallow evo¬ 
lution of f e sc£ion and f e action vs as a function of dust 
extinction for the density-bounded nebula model (Model 
B). Over the dynamic range in /3, the ModelB_dust model 
evolves « 0.5dex from log 10 / eS c£i°n(,d = —2.3) « 24.75 to 
log 10 fe S c&on{/3 = -1.7) RS 24.25. 

The second relation, ModelAJesc, follows the / esc £ion or 
/esc«ion vs /3 evolution as a function of / esc for the ionization 
bounded nebula model (Model A). This model evolves from 
logic /esc£ion « 25 at /? = —2.3 to effectively zero ionizing 
photons per unit luminosity/star-formation at /3 = —1.8. 
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Due to the lack of constraints on ft for galaxies all the way 
down to Muv = —13, we set a lower limit on how steep 
the UV slope can become. This limit is chosen to match 
the UV slope for the dust-free, / esc = 1 scenario for fiducial 
model and has a slope of ft = —2.55. We choose these two 
relations (three including the constant assumption above) 
because they correspond to the most likely mechanisms 
through which ft or the LyC escape fraction are expected 
to evolve. 

Firstly, the constant (equivalent to ‘ModelA_dust’) and 
‘ModelB-dust’ models cover the assumption that evolution 
in the dust content of galaxies is responsible for the ob¬ 
served evolution in ft and any corresponding evolution in 
the ionizing efficiency of galaxies (/ e scfion or / es cKion). Sec¬ 
ondly, the ‘ModelA_fesc ’ model corresponds to an evolution 
in /esc alone and for the ionization-bounded nebula with 
holes model, represents the steepest evolution of / eS c£ion 
(//escftion) with respect to ft of any of the parameters. While 
there is no obvious physical mechanism for such evolution 
at high-redshift, using this model we can at least link an in¬ 
ferred /esc redshift evolution such as that from the Kuhlen & 
Faucher-Giguere (2012) and R13 to a corresponding evolu¬ 
tion in ft and vice-versa. In these works, the evolving escape 
fraction is parametrised as: 

(14) 

where fo = 0.054 and 7 = 2.4 (R13) and are constrained 
by the observed IGM emissivity values of Faucher-Giguere 
et al. (2008) at 2 ^ 4 and the WMAP total integrated opti¬ 
cal depth measurements of Hinshaw et al. (2013) at higher 
redshift. For the subsequent analysis we also show how the 
total ionizing emissivity would change following this evolu¬ 
tion of / esc (assuming a constant £i on = 25.2 in line with 
R13). 

It is important to note here that the more recent mea¬ 
surements of the IGM emissivity at 2 ~ 4 by Becker & 
Bolton (2013) are a factor ~ 2 greater than those of Faucher- 
Giguere et al. (2008) at the same redshift. As such, the 
assumed values may under-estimate the zero-point fo and 
over-estimate the slope of the / eac evolution compared to 
those fitted to the IGM emissivities of Becker & Bolton 
(2013). However, we include this to illustrate the effects 
that forcing consistency with IGM and optical depth mea¬ 
surements has on total ionizing emissivity for comparable 
underlying UV luminosity/star-formation rate density mea¬ 
surements relative to the assumption of a constant conver¬ 
sion. 

For the star-formation history assumed throughout this 
work, both of these models should in principle also take into 
account the reddening of the intrinsic UV slope due to age 
evolution (see Section 3.4). However, we neglect this contri¬ 
bution in the following analysis for two reasons. Firstly, the 
effect of age evolution on the observed ft is small in com¬ 
parison to that of /esc or dust for these models. Secondly, 
the vector relating ft and f e sc £ ion f° r increasing age is either 
almost parallel to or slightly steeper than those for / esc and 
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Figure 9. Ionizing emissivity, Ni ou , predicted by the luminos¬ 
ity functions measured by Bouwens et al. (2014a) for an evolving 
/escCion as a function of redshift, based on the luminosity weighted 
average ft . See text for details on the assumed / es c£ion as a, func¬ 
tion of (/3) uv f° r the ModelB_dust (top - green) and ModelA_fesc 
(bottom - blue) models. In both panels, the thick solid and dashed 
lines correspond to the UV luminosity density integrated down to 
the observational limit and a constant Muv = —13 respectively. 
The yellow lines indicate the emissivity predicted for redshift de¬ 
pendent /esc as inferred by R13, see Equation 14. As in Figure 8, 
we show the IGM emissivity measurements and corresponding to¬ 
tal errors of Becker & Bolton (2013) (tan line and cross-hatched 
region respectively) and the ionizing background constraints of 
Bouwens et al. (2015) (blue-gray diagonal-hatched region). 

dust respectively. As such, the effects on the inferred / eS c£ ion 
will be negligible. 

ft.3.1 Evolving luminosity-averaged f esc £ion 

To explore how a ft dependent coefficient would change the 
inferred emissivities, we first calculate a constant / eac £ion for 
each redshift based on the luminosity weighted {ft) vv . The 
2 ~ 4 to z ~ 7 (/3} uv used are those from Bouwens et al. 
(2014b), with the 2 ~ 8 value based on the fit outlined in 
Equation 7. 

For clarity, we plot the resulting predicted emissivities 
only for the UV luminosity densities predicted by Bouwens 
et al. (2014a) luminosity function parametrisations, these 
are shown in Fig. 9 for both the ModelB_dust and Mod- 
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elA_fesc (top and bottom panels respectively) / esc £ion 03) 
relations. The shaded regions around the solid green and 
blue lines (top and bottom panels respectively) represent 
the uncertainty on / esc £ion due to the statistical uncertainty 
in (/3) vv . The full statistical uncertainties include the lumi¬ 
nosity density errors illustrated in Fig. 8 and are included 
in Appendix Table 3. 

By assuming a /3 dependent coefficient, the estimated 
Nion evolution changes shape to a much shallower evolution 
with redshift. For our ModelB_dust j3 relation, the decline 
in JVion from 2 = 4 to 2 = 8 is reduced by tts 0.25 dex when 
the luminosity function is integrated down to the limit of 
Muv < —13. The effect is even stronger for the ModelA_fesc 
/3 evolution, with the Mon actually increasing over this red- 
shift. The larger / eac £ion inferred by the increasingly blue UV 
slopes at high redshift are able to balance the rapid decrease 
in UV luminosity density. 

A key effect of the increasing ionizing efficiency with 
increasing redshifts is that for both f esc ^ion(/3) relations ex¬ 
plored here, the observable galaxy population at 2 ~ 7 is 
now capable of maintaining reionization for a clumping fac¬ 
tor of Chu = 3. This is in contrast to the result inferred 
when assuming a constant f e sc £km- 

However, we know the brightest galaxies are in fact also 
the reddest and that any of our predicted f e sc £ion(/3) rela¬ 
tions imply they are therefore the least efficient at produc¬ 
ing ionizing photons. The application of an average / eS c£ion 
(even one weighted by the relative contributions to the lu¬ 
minosity density) may give a misleading impression of the 
relative contribution the brightest galaxies make to the ion¬ 
izing background during the EoR. A more accurate picture 
can be obtained by applying a luminosity dependent /3 re¬ 
lation (e.g. /esc£ion(A/uv)) to observed luminosity function 
and integrating the ionizing emissivity, Ni on , from this. 


4-3.2 Luminosity dependent f eac £ion 

Using the observed /?(Muv) relations of Bouwens et al. 
(2014b) (Fig. 3) and our models for / eS c£ion(/3), we next 
calculate Mj on a function of both the changing luminos¬ 
ity function and the evolving colour magnitude relation. 
Fig. 10 shows the evolution of Mon based on these assump¬ 
tions, again for the Bouwens et al. (2014a) luminosity func¬ 
tion parametrisations. For both the ModelB-dust and Mod- 
elA_fesc relations, the emissivity of galaxies above the lim¬ 
iting depths of the observations are reduced due to the fact 
that the brighter galaxies have significantly redder observed 
/3s. 

I 11 the case of the ModelA_fesc /3 evolution and the 
ModelB_dust relation at high redshift, the difference be¬ 
tween the Mon for Muv < —17 and Muv < —13 is quite 
significant. This is due to the observed steepening of the 
colour-magnitude relation (Section 2.3) at higher redshifts 
results in a larger / esc £i on evolution between the brightest 
and faintest galaxies in the luminosity function. When inte¬ 
grating the UV luminosity function from fainter magnitudes, 
the number of ionizing photons produced per unit UV lumi- 
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Figure 10. Ionizing emissivity, 2Vi on , predicted by the luminos¬ 
ity functions measured by Bouwens et al. (2014a) for a lumi¬ 
nosity dependent / es c£ion- See text for details on the assumed 
fe sc£ion as a function of /3(Muv) f° r the ModelB_dust( top) and 
ModelA_fesc (bottom) models. In both panels, the thick solid, 
dashed and dotted lines correspond to the UV luminosity density 
integrated down to Muv = —17, -15 and -13 respectively. The 
yellow lines indicate the emissivity predicted for redshift depen¬ 
dent /e S c as inferred by R13, see Equation 14. As in Figure 8, we 
show the IGM emissivity measurements and corresponding to¬ 
tal errors of Becker & Bolton (2013) (tan line and cross-hatched 
region respectively) and the ionizing background constraints of 
Bouwens et al. (2015) (blue-gray diagonal-hatched region). 


nosity density significantly increases. At lower redshifts, the 
sharp decline in / esc £ion for redder galaxies in ModelAJesc 
means that the brightest galaxies can potentially contribute 
very little to the total ionizing emissivity. 

The total galaxy ionizing emissivity (Muv < —13) for 
both models of /3 evolution is high enough to maintain reion¬ 
ization (for a clumping factor of Chii < 3) at all redshifts. 
In fact, only galaxies down to a rest-frame magnitude of 
Muv = —15 are required to match these rates. Furthermore, 
despite the increased Mon predicted at 2 > 5 both models 
are in good agreement with the observed IGM emissivities of 
Becker & Bolton (2013) at lower redshifts when the luminos¬ 
ity functions are integrated down to a limit of Muv < —13. 
Including galaxies such faint galaxies does however result 
in a potential over-production of ionizing photons at 2 > 6 
based on the emissivity estimates of Bouwens et al. (2015). 
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Given that we have shown that significantly lower values of 
/esc are still consistent with the observed colours (as stated 
in Section 3.2), an overproduction for the optimistic models 
assumed does not necessarily present an immediate problem. 

From these results we can see that changes in the ion¬ 
izing efficiency of galaxies during EoR which are still con¬ 
sistent with the evolving UV continuum slopes have signif¬ 
icantly less effect on the predicted total ionizing emissivity 
(Muv < —13) at z ~ 4 than at higher redshifts based on the 
current observations. This is a crucial outcome with regards 
to current numerical models for the epoch of reionization, as 
it allows for a wider range of reionization histories which are 
still consistent with both the observed UV luminosity/SFR- 
density and lower redshift IGM emissivity estimates. 


5 DISCUSSION AND FUTURE PROSPECTS 

In several previous studies of the reionization history of 
the Universe the conclusion has been drawn that at earlier 
times in the epoch of reionization, galaxies must have been 
more efficient at ionizing the surrounding IGM than similar 
galaxies at lower redshift (Becker & Bolton 2013; Kuhlen & 
Faucher-Giguere 2012; Robertson et al. 2013). Based on the 
constraints on galaxy stellar populations and escape frac¬ 
tions explored in Section 3 and their application to the ex¬ 
isting observations in Section 4, it is not yet possible to 
establish that one particular galaxy property is evolving to 
cause such an increase in ionizing efficiency. 

However, what we find in this work is that evolution 
in galaxy properties, such as dust extinction and escape 
fraction (or some combination of these and others), which 
are consistent with the observed colour evolution of high- 
redshift galaxies can readily account for any increase in 
the ionizing efficiency required by other constraints such as 
the total optical depth. Ongoing and future observations of 
both local and distant galaxies will be able to provide much 
tighter constraints on the evolving galaxy properties. 

If the observed /3 evolution is a result of dust alone, as 
is assumed by Bouwens et al. (2012b) and other works, the 
inferred evolution in / es cCion as a function of ft is highly de¬ 
pendent on the assumed model of Lyman continuum escape 
and hence the underlying geometry of dust and gas. For ex¬ 
ample, if the channels through which LyC photons escape 
are dust-free (as in Model A), the effect of the dust evolu¬ 
tion will have negligible effect on the emissivity of galaxies 
as a function of fi (as discussed in Section 3.4). Real galaxies 
will of course be significantly more complicated (and messy) 
than the simple toy models adopted in this work, as such the 
channels through which LyC escape may also contain signif¬ 
icant quantities of dust. We find that if we modify Model A 
such that the dust screen is extended to include the low- 
density channels, the resulting model is indistinguishable 
from Model B (density bounded nebula) with regards to 
/3 as a function of / esc or Av . Such a model closely matches 
that observed by Borthakur et al. (2014) (see also Heckman 
et al. 2011) for a local analog of galaxies during the EoR 
and represents our best model for LyC escape. 


For any of the plausible causes for the luminosity and 
redshift evolution of ft (dust, metallicity, escape fraction), 
the models explored in this work infer that fainter/low- 
mass galaxies are emitting more ionizing photons per unit 
star-formation into the IGM than their higher mass coun¬ 
terparts. Currently, simulations of galaxies at high redshift 
draw somewhat differing conclusions on the mass/luminosity 
dependence of the escape fraction. Based on a combination 
of theoretical models and the existing limited observations, 
Gnedin et al. (2008) find that angular averaged escape frac¬ 
tion increases with higher star-formation rates and galaxy 
masses, the inverse of what we predict based on j3 alone. 
However, in isolation from the model predictions, the ob¬ 
servational data explored by Gnedin et al. (2008) does not 
place any strong constaints on the luminosity dependence 
of /esc (Giallongo et al. 2002; Fernandez-Soto et al. 2003; 
Shapley et al. 2006). 

Subsequent simulations predict the opposite luminosity 
dependence, in better agreement with the colour evolution 
predictions of this work (Razoumov & Sommer-Larsen 2010; 
Yajima et al. 2010). Recent work exploring the escape frac¬ 
tion of both typical (Kimm & Cen 2014) and dwarf (Wise 
et al. 2014) galaxies at z ^ 7 find that the instantaneous 
escape fraction (measured at the virial radius in these simu¬ 
lations) is inversely proportional to the halo mass. However, 
as discussed by Kimm & Cen (2014), the average instan¬ 
taneous escape fraction may be somewhat misleading due 
to the bursty nature of star-formation in their models and 
the delay between the peak SFR and maximum escape frac¬ 
tion for an episode of star-formation. They find that the 
time-averaged escape fraction weighted by the overall LyC 
photon production rate remains roughly constant for size 
haloes. Improved measurements on the stellar or halo mass 
dependence of / esc are therefore clearly crucial. 

Although direct measurements of the LyC escape frac¬ 
tion for galaxies during EoR will never be possible due to 
the effects of IGM absorption along the line of sight, mea¬ 
surements of the escape fraction as a function of stellar mass 
and luminosity (or SFR) at z < 3 should soon be possible 
due to the deep UV imaging of new surveys such as the 
UVUDF (Teplitz et al. 2013) and the forthcoming GOODS 
UV Legacy Survey (PI: Oesch, G013872). The wealth of an¬ 
cillary data available in these fields (both photometric and 
spectroscopic) should make it possible to tightly constrain 
/esc, Cion or Ki on , and /3 for either individual galaxies or sam¬ 
ples stacked by galaxy properties. Measuring /3 vs / esc £ion 
at 2 < z < 3 would significantly reduce systematic errors in 
the inferred Ni on during the EoR from incorrect or poorly 
informed assumptions on / ea cCion- 

Given the large degeneracies in /? with respect to the 
various stellar population parameters, understanding /3 vs 
/escCion both at z ~ 3 and during the EoR will require an 
improved understanding of what is responsible for the ob¬ 
served P evolution. With ALMA observations of statistically 
significant samples of galaxies at high-redshift, it should be 
possible to make strong constraints on not just the strength 
and attenuation curve of the dust extinction, but also the 


© 2015 RAS, MNRAS 000, i-22 


20 Duncan and Conselice 


location and geometry of the dust relative to the gas and 
star-formation within galaxies (De Breuck et al. 2014). 

With the new generation of near-infrared sensitive spec¬ 
trographs allowing precision spectroscopic measurement of 
metallicities and dust out to 2 > 3 (e.g. MOSFIRE: Kriek 
et al. (2014)), it will be possible to place much more ac¬ 
curate priors on the expected ages, metallicities and star- 
formation histories for galaxies during the EoR. Finally, as 
with many outstanding problems in astrophysics, the launch 
of the James Webb Space Telescope will address many of the 
systematic and statistical uncertainties which limit current 
observations. Crucially, JWST should be able to probe much 
fainter galaxy populations, potentially down to rest-frame 
magnitudes of Afuv = —15 and below. Based on the find¬ 
ings in this paper, such observations might even mean we are 
finally able to observe the full galaxy population responsible 
for powering reionization. 


6 SUMMARY 

In this work, we explore in-depth the ionizing photon budget 
of galaxies during the epoch of reionization based solely on 
the observed galaxy properties. For the latest observational 
constraints on the star-formation rate and UV luminosity 
density at 2 > 4, we assess the ionizing emissivity consis¬ 
tent with new constraints on the rest-frame UV colours of 
galaxies at these redshifts. 

Using a comprehensive set of SED models for two plau¬ 
sible Lyman continuum escape mechanisms - previously out¬ 
lined in Zackrisson et al. (2013) - we explore in detail the re¬ 
lationship between the UV continuum slope /5 and the num¬ 
ber of ionizing photons produced per unit UV luminosity 
or star-formation (£i on and Ki on respectively). We find that 
the ionizing efficiencies assumed by several previous works 
(login / eS c£ion = 24.5 — 24.6: Robertson et al. (2013); Kuhlen 
& Faucher-Giguere (2012)) are still consistent with the cur¬ 
rent ft observations during the EoR. However, for both of 
the LyC escape models explored here, this assumption is 
close to the maximum efficiency which is still consistent with 
the fiducial UV slope typically considered at these redshifts 
(/3 = —2). Based on our SED modelling, escape fractions or 
ionizing efficiency which are 1 dex lower than typically as¬ 
sumed are still consistent with the observed galaxy colours. 

Applying the fiducial log 10 / esc £ion = 24.5 to the the 
latest observations of the luminosity and mass functions at 
2^4, we find that at 2 ~ 6, the observed population can 
produce enough ionizing photons to maintain reionization 
assuming a clumping factor Chu = 3. At earlier times, we 
confirm previous results which found that galaxies from be¬ 
low our current observation limits are required to produce 
enough ionizing photons to maintain reionization at 2 ~ 7 
and beyond (Robertson et al. 2010, 2013; Kuhlen & Faucher- 
Giguere 2012; Finkelstein et al. 2012b). 

Motivated by increasing evidence for a luminosity de¬ 
pendence of the UV continuum slope (Bouwens et al. 2014b; 
Rogers et al. 2014) and evidence for evolution in this relation 


with redshift (Bouwens et al. 2014b), we explore the effects 
of assuming an ionizing efficiency which is not constant but 
varies with the observed /3. The two galaxy properties that 
are able to plausibly account for the required range of ob¬ 
served /3’s, dust extinction and /esc, both predict an ionizing 
efficiency which increases for increasingly blue UV contin¬ 
uum slopes. While the other galaxy properties such as age 
and metallicity predict similar trends, current observations 
do not support a large enough variation to account for the 
required range of observed UV slopes. 

We find that when assuming an ionizing efficiency based 
on the luminosity-weighted average /3, the currently observ¬ 
able galaxy population alone is now able to maintain reion¬ 
ization at 2 ~ 7 (assuming CW = 3). Despite this increase in 
efficiency at early times, the predicted Ni on at 2 < 5 remain 
consistent with measurements based on the IGM (Becker & 
Bolton 2013). 

Assuming instead that the ionizing efficiency of galaxies 
is dependent on their luminosity, the observed Muv — P re¬ 
lations and our SED models can result in significant changes 
in the inferred ionizing photon budget. Since our models sug¬ 
gest that redder (brighter) galaxies have lower / es c£km than 
their blue (faint) counterparts, the inferred ionizing photon 
budget for the currently observable galaxy population may 
be significantly reduced, especially at lower redshifts. How¬ 
ever, because of the increasing importance of faint galaxies 
(which have higher inferred / eS c£ion), only galaxies down to 
Muv ~ —15 may be required to produce the required ioniz¬ 
ing photons. 

Our conclusion is that the inferred ability of galaxies 
to complete or maintain reionization is highly dependent on 
the stellar population assumptions used to predict their ion¬ 
izing efficiencies. Crucially though, the models explored in 
this study can potentially allow for a wide range of reioniza¬ 
tion histories whilst remaining consistent with the observed 
colour evolution and luminosity (or star-format ion rate) den¬ 
sity during this epoch. Future work on constraining both the 
colour and luminosity dependence of / esc at lower redshifts 
as well as measuring the ages and dust content of galaxies 
during the EoR will be vital in understanding the precise 
ionizing emissivity of galaxies throughout this epoch. 
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APPENDIX A: DATA TABLES 

In this appendix we present the estimated puv and Mon at 
4 ^ 2 : ^ 8 for the range of assumptions outlined in Section 4. we 
list the calculated properties for integration limits M\ jy = -17, -15 
and -13 respectively. In Table A we list the values based on the lu¬ 
minosity function parametrisations of Bouwens et al. (2014a) and 
are plotted in Figures 8, 9 and 10. In Table A, we list the corre¬ 
sponding values for the luminosity function parametrisations of 
(Finkelstein et al. 2014). For both sets of values, we include errors 
based on the uncertainties in the luminosity function parameters 
and the random errors in the weighted average of /3 or the best-fit 
(3 — Muv slope parameters (Bouwens et al. 2014b) as appropri¬ 
ate. When assuming an SMC-like dust attenuation (Pei 1992), 
to match the fiducial (3 = — 2 and log 10 fe sc£ion = 24.5 for LyC 
escape Model B, we find f esc =0.28 and Ay =0.08 are required. 
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Table Al. Calculated values of puv and Mon f° r the different integration limits and efficiency assumptions explored in the paper, based 
on the luminosity function parametrisations of Bouwens et al. (2014a). For each calculated value, we include statistical errors from the 
uncertainties in the Schechter (1976) parameters and (3 observations. Also shown are the effects of some of the assumptions made in 
Section 3.1 and their corresponding systematic changes to the estimated values. 
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Table A2. Calculated values of puv and Mon f° r the different integration limits and efficiency assumptions explored in the paper, based 
on the luminosity function parametrisations of Finkelstein et al. (2014). For each calculated value, we include statistical errors from the 
uncertainties in the Schechter (1976) parameters and (3 observations. Also shown are the effects of some of the assumptions made in 
Section 3.1 and their corresponding systematic changes to the estimated values. 
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